###############################
# This script produces Figure 2 in the Paper -- Monte Carlo Results
# It relies on the "Functions_Mode_Agree_Boot.R" script and draws them in from same directory
# Libraries are loaded in script 1. 
# Last Updated 6 September 2018

rm(list=ls())
set.seed(1)


source('Functions_Mode_Agree_Boot.R')

#################
### Create a scenario with bimodal parties --- two different versions of reality
### For a subset of extreme parties, experts sometimes view mirror image of the party

partyrev <- function (partypos){

	extreme <- rep(NA, length(partypos))

	extreme[partypos>7.5|partypos<2.5] <- 1

	index <- c(1:length(partypos))

	samp <- sample(index[is.na(extreme)!=1] , round(0.35*sum(extreme,na.rm=T),0))

	partypos.rev <- partypos

	partypos.rev[samp] <- 10-partypos[samp]

	return(partypos.rev)

}



### Create a function to generate expert assessments

ex.assess <- function(nparty,nexpert,partypos,noise,stretch,shift,rev=FALSE){
		
	# Create holding matrix
	exdta <- matrix(NA,nrow=nparty,ncol=nexpert)
	
	# generate random expert-specific bias
	stretch.bias <- 1+runif(nexpert,-1*stretch,stretch)
	shift.bias <- runif(nexpert,-1*shift,shift)
	
	# generate party-specific standard dev for noise paramenter
	sd <- runif(nparty,0,noise)
	
	if (rev==TRUE) { 
		
		partypos.rev <- partyrev(partypos) 
		
		for (i in 1:nparty){
			for (j in 1:nexpert){
			
				if (runif(1) > 0.65 ) {
				
					exdta[i,j] <- partypos.rev[i]*stretch.bias[j]+shift.bias[j]+rnorm(1,0,sd[i])	

				}
			
			
			else{
			
					exdta[i,j] <- partypos[i]*stretch.bias[j]+shift.bias[j]+rnorm(1,0,sd[i])	
			
				}
			}
		}
	}
	
	else{
		
		for (i in 1:nparty){
	
			exdta[i,] <- partypos[i]*stretch.bias+shift.bias+rnorm(nexpert,0,sd[i])	
	
		}
	
	}
		
		
	exdta <- round(exdta,0)

	exdta[exdta > 10] <- 10
	exdta[exdta < 0] <- 0
	
	return(exdta)	

}


#################################

mad.lm <- function (expert.dta) {

	mad.mean <-NULL
	mad.median <-NULL
	mad.mode <-NULL

	coef.mean <- NULL
	coef.median <- NULL
	coef.mode <- NULL

	for (i in 1: 1000){
	
		mexp <- apply(expert.dta[[i]],1,mean)
		mexp.1 <- apply(expert.dta[[i]],1,median)
		mexp.2 <- apply(expert.dta[[i]],1,Mode)
		
		mad.mean[i] <- mean(abs(mexp-partypos))
		mad.median[i] <- mean(abs(mexp.1-partypos))
		mad.mode[i] <- mean(abs(mexp.2-partypos))
		
		coef.mean[i] <- lm(mexp~partypos)[[1]][2]
		coef.median[i] <- lm(mexp.1~partypos)[[1]][2]
		coef.mode[i] <- lm(mexp.2~partypos)[[1]][2]
	
	}
	
	out.mat <- cbind(mad.mean, mad.median, mad.mode, coef.mean, coef.median, coef.mode)
	return(out.mat)
}


################################################
#### Monte Carlo Comparing Bias of Mean, Median and Mode
################################################


### Latent Party Positions -- 100 parties ranging from 0 to 10, 15 experts
nparty <- 100
partypos <- runif(nparty,0,10)

### Set noise stretch and shift
noise <- 3
stretch <- 0.3
shift <- 2



###################################
#### Run Monte Carlo -- 1000 times -- with different numbers of experts

expert.dta.5 <- vector("list",length=1000)
expert.dta.10 <- vector("list",length=1000)
expert.dta.15 <- vector("list",length=1000)
expert.dta.20 <- vector("list",length=1000)

for (i in 1:1000) {
	expert.dta.5[[i]] <- ex.assess(nparty,5,partypos,noise,stretch,shift)
	expert.dta.10[[i]] <- ex.assess(nparty,10,partypos,noise,stretch,shift)
	expert.dta.15[[i]] <- ex.assess(nparty,15,partypos,noise,stretch,shift)
	expert.dta.20[[i]] <- ex.assess(nparty,20,partypos,noise,stretch,shift)
}


outexp.5  <- mad.lm(expert.dta.5)
outexp.10 <- mad.lm(expert.dta.10)
outexp.15 <- mad.lm(expert.dta.15)
outexp.20 <- mad.lm(expert.dta.20)

####

expert.dta.5.bi <- vector("list",length=1000)
expert.dta.10.bi <- vector("list",length=1000)
expert.dta.15.bi <- vector("list",length=1000)
expert.dta.20.bi <- vector("list",length=1000)

for (i in 1:1000) {
	expert.dta.5.bi[[i]] <- ex.assess(nparty,5,partypos,noise,stretch,shift,rev=T)
	expert.dta.10.bi[[i]] <- ex.assess(nparty,10,partypos,noise,stretch,shift,rev=T)
	expert.dta.15.bi[[i]] <- ex.assess(nparty,15,partypos,noise,stretch,shift,rev=T)
	expert.dta.20.bi[[i]] <- ex.assess(nparty,20,partypos,noise,stretch,shift,rev=T)
}

outexp.5.bi  <- mad.lm(expert.dta.5.bi)
outexp.10.bi <- mad.lm(expert.dta.10.bi)
outexp.15.bi <- mad.lm(expert.dta.15.bi)
outexp.20.bi <- mad.lm(expert.dta.20.bi)

### col order from output: mad.mean, mad.median, mad.mode, coef.mean, coef.median, coef.mode

coef.mean <- cbind(c(rep(5,1000),rep(10,1000),rep(15,1000),rep(20,1000)),c(outexp.5[,4],outexp.10[,4],outexp.15[,4],outexp.20[,4]))
coef.median <- cbind(c(rep(5,1000),rep(10,1000),rep(15,1000),rep(20,1000)),c(outexp.5[,5],outexp.10[,5],outexp.15[,5],outexp.20[,5]))
coef.mode <- cbind(c(rep(5,1000),rep(10,1000),rep(15,1000),rep(20,1000)),c(outexp.5[,6],outexp.10[,6],outexp.15[,6],outexp.20[,6]))

mad.mean <- cbind(c(rep(5,1000),rep(10,1000),rep(15,1000),rep(20,1000)),c(outexp.5[,1],outexp.10[,1],outexp.15[,1],outexp.20[,1]))
mad.median <- cbind(c(rep(5,1000),rep(10,1000),rep(15,1000),rep(20,1000)),c(outexp.5[,2],outexp.10[,2],outexp.15[,2],outexp.20[,2]))
mad.mode <- cbind(c(rep(5,1000),rep(10,1000),rep(15,1000),rep(20,1000)),c(outexp.5[,3],outexp.10[,3],outexp.15[,3],outexp.20[,3]))


coef.mean.bi <- cbind(c(rep(5,1000),rep(10,1000),rep(15,1000),rep(20,1000)),c(outexp.5.bi[,4],outexp.10.bi[,4],outexp.15.bi[,4],outexp.20.bi[,4]))
coef.median.bi <- cbind(c(rep(5,1000),rep(10,1000),rep(15,1000),rep(20,1000)),c(outexp.5.bi[,5],outexp.10.bi[,5],outexp.15.bi[,5],outexp.20.bi[,5]))
coef.mode.bi <- cbind(c(rep(5,1000),rep(10,1000),rep(15,1000),rep(20,1000)),c(outexp.5.bi[,6],outexp.10.bi[,6],outexp.15.bi[,6],outexp.20.bi[,6]))

mad.mean.bi <- cbind(c(rep(5,1000),rep(10,1000),rep(15,1000),rep(20,1000)),c(outexp.5.bi[,1],outexp.10.bi[,1],outexp.15.bi[,1],outexp.20.bi[,1]))
mad.median.bi <- cbind(c(rep(5,1000),rep(10,1000),rep(15,1000),rep(20,1000)),c(outexp.5.bi[,2],outexp.10.bi[,2],outexp.15.bi[,2],outexp.20.bi[,2]))
mad.mode.bi <- cbind(c(rep(5,1000),rep(10,1000),rep(15,1000),rep(20,1000)),c(outexp.5.bi[,3],outexp.10.bi[,3],outexp.15.bi[,3],outexp.20.bi[,3]))


# Use pdf to save plot if required.
pdf("MCHistNew.pdf",width=10,height=7)

par(mfrow=c(2,3),oma=c(0,11,0,0))

boxplot(coef.mean[,2]~coef.mean[,1], ylim = c(0.4,1.2),ylab="coefficients", xlab="number of experts", main="Mean")
abline(h=1, lty=3)
boxplot(coef.median[,2]~coef.median[,1], ylim = c(0.4,1.2),xlab="number of experts", main="Median")
abline(h=1, lty=3)
boxplot(coef.mode[,2]~coef.mode[,1], ylim = c(0.4,1.2),xlab="number of experts", main="Mode")
abline(h=1, lty=3)

boxplot(coef.mean.bi[,2]~coef.mean.bi[,1],ylim = c(0.3,1.15), ylab="coefficients",xlab="number of experts", main="")
abline(h=1, lty=3)

boxplot(coef.median.bi[,2]~coef.median.bi[,1],ylim = c(0.3,1.15), xlab="number of experts",main="")
abline(h=1, lty=3)

boxplot(coef.mode.bi[,2]~coef.mode.bi[,1],ylim = c(0.3,1.15),xlab="number of experts", main="")
abline(h=1, lty=3)

mtext("Unimodal",side=2,outer=T,las=1,adj=0,line=10,at=0.75)
mtext("Bimodal",side=2,outer=T,las=1,adj=0,line=10,at=0.25)

dev.off()